Lab sessions: Day 1
Rob J Hyndman
25 June 2018
Lab Session 1
retaildata <- read.csv("retail.csv")
mytimeseries <- ts(retaildata[,4], frequency=12, start=c(1982,4))
autoplot(mytimeseries)
Lab Session 2
There is a strong trend to around 2008, weak seasonality (with a jump in December each year), and some evidence of cycles (with dips in 2001 and 2012).
Lab Session 3
bicoal
bicoal is annual data, so you can’t do the seasonal plots.
dole
Other series are similar.
Lab Session 4
##
## Ljung-Box test
##
## data: Residuals from Seasonal naive method
## Q* = 32.269, df = 8, p-value = 8.336e-05
##
## Model df: 0. Total lags used: 8There is some remaining autocorrelation in the residuals: the Null of no joint autocorrelation is clearly rejected. We can also see a significant spike on the seasonal (4th lag) in the ACF. There is considerable information remaining in the residuals which has not been captured with the seasonal naïve method. The residuals do not appear to be too far from Normally distributed.
Lab Session 5
train <- window(mytimeseries, end=c(2010,12))
test <- window(mytimeseries, start=2011)
autoplot(cbind(Training=train,Test=test))
## ME RMSE MAE MPE MAPE MASE
## Training set 5.502402 27.74935 19.33784 3.369450 10.447161 1.0000000
## Test set 5.227586 23.40458 18.60000 1.619239 8.172045 0.9618449
## ACF1 Theil's U
## Training set 0.8703252 NA
## Test set 0.4544630 0.9518744The number to look at here is the test set RMSE. That provides a benchmark for comparison when we try other models.
##
## Ljung-Box test
##
## data: Residuals from Seasonal naive method
## Q* = 1256.8, df = 24, p-value < 2.2e-16
##
## Model df: 0. Total lags used: 24The residuals do not look like white noise there are lots of dynamics left over that need to be explored. They also do not look close to normal, with very long tails.
The accuracy measure are always sensitive to the training/test split. There are better ways to check the robustness of the methods in terms of accuracy such as using a tsCV().
Lab Session 6
e_mean <- tsCV(mytimeseries, meanf)
e_naive <- tsCV(mytimeseries, naive)
e_snaive <- tsCV(mytimeseries, snaive)
e_drift <- tsCV(mytimeseries, rwf, drift=TRUE)
# Construct squared CV errors matrix
e2 <- cbind(Mean=e_mean^2, Naive=e_naive^2,
SNaive=e_snaive^2, Drift=e_drift^2)
# Remove rows with any missing for a fair comparison
e2 <- na.omit(e2)
# Find MSE values
colMeans(e2)## Mean Naive SNaive Drift
## 4685.8665 410.0438 691.7144 411.4290Lab Session 7
## ME RMSE MAE MPE MAPE MASE
## Training set -2.740622 26.56395 19.38206 -2.872958 10.05648 0.9560879
## ACF1
## Training set -0.005983602## ME RMSE MAE MPE MAPE MASE
## Training set 0.04499087 26.58219 19.18491 -1.142201 9.653791 0.9463626
## ACF1
## Training set 0.01348202## ME RMSE MAE MPE MAPE MASE
## Training set -2.891496 26.54019 19.2795 -2.907633 10.01894 0.9510287
## ACF1
## Training set -0.003195358These RMSE values are not really comparable because the models have different numbers of parameters.
The best model is the last one.
##
## Ljung-Box test
##
## data: Residuals from Damped Holt's method
## Q* = 10.782, df = 5, p-value = 0.05587
##
## Model df: 5. Total lags used: 10The residuals are pretty good and (just) pass the Ljung-Box test. However, the forecast make little sense (with the prediction intervals going negative very quickly).
Lab Session 8
The seasonal variation increases with the level of the series. So we need to use multiplicative seasonality.
## ME RMSE MAE MPE MAPE MASE
## Training set -0.2054187 9.280748 6.244629 -0.2790438 3.338901 0.3336358
## ACF1
## Training set 0.05182697## ME RMSE MAE MPE MAPE MASE
## Training set 0.5050393 9.194427 6.128123 0.2458937 3.290157 0.3274111
## ACF1
## Training set -0.04640084There is not much difference between these models, and we would expect the trend to continue, so I would prefer to use the non-damped version.
##
## Ljung-Box test
##
## data: Residuals from Holt-Winters' multiplicative method
## Q* = 41.483, df = 8, p-value = 1.693e-06
##
## Model df: 16. Total lags used: 24There are significant correlations in the residuals, especially at lag 12. So these residuals do not look like white noise.
train %>% window(end=c(2010,12)) %>%
hw(seasonal='multiplicative', damped=FALSE) %>%
accuracy(x=mytimeseries)## ME RMSE MAE MPE MAPE MASE
## Training set -0.1212843 9.260807 6.054543 -0.03697659 3.442954 0.3130931
## Test set 10.1364618 19.695875 15.215508 4.30118709 6.977221 0.7868257
## ACF1 Theil's U
## Training set 0.09051818 NA
## Test set 0.61161198 0.909366The test set RMSE is much better than for the seasonal naïve method.
Lab Session 9
## ETS(M,Ad,M)
##
## Call:
## ets(y = mytimeseries)
##
## Smoothing parameters:
## alpha = 0.7718
## beta = 0.0025
## gamma = 0.0561
## phi = 0.98
##
## Initial states:
## l = 64.0899
## b = 0.2813
## s = 0.9869 0.9362 1.0043 1.1705 1.0216 1.0293
## 0.9848 0.9998 0.993 0.9436 0.9694 0.9606
##
## sigma: 0.0464
##
## AIC AICc BIC
## 4404.041 4405.697 4477.272
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.5339249 9.613425 6.241919 0.1970991 3.257992 0.333491
## ACF1
## Training set -0.03180379This is equivalent to a damped-trend multiplicative Holt-Winters’ method. The very small \(\beta\) and \(\gamma\) values show that the seasonality and trend change slowly. There is also very little damping (\(\phi\) is close to 1).
Not too bad given the noisy data.
These are great forecasts which have adapted to the trend and seasonality really well.
The seasonality seems to be over-done, but the wide prediction intervals are probably sensible apart from the fact that they go negative.
These data are cyclic, and ETS does not handle cycles in time series.
Lab Session 10
train <- window(mytimeseries, end=c(2010,12))
test <- window(mytimeseries, start=2011)
f1 <- snaive(train, h=length(test))
f2 <- rwf(train, h=length(test))
f3 <- rwf(train, drift=TRUE, h=length(test))
f4 <- meanf(train, h=length(test))
f5 <- hw(train, h=length(test),
seasonal='multiplicative')
f6 <- ets(train) %>% forecast(h=length(test))
c(
SNaive=accuracy(f1, test)[2,"RMSE"],
Naive=accuracy(f2, test)[2,"RMSE"],
Drift=accuracy(f3, test)[2,"RMSE"],
Mean=accuracy(f4, test)[2,"RMSE"],
HW=accuracy(f5, test)[2,"RMSE"],
ETS=accuracy(f6, test)[2,"RMSE"])## SNaive Naive Drift Mean HW ETS
## 23.40458 37.09060 56.64803 58.73867 59.85277 25.43061These results might depend on the peculiarities of the test set. That’s why we usually prefer to do a cross-validation approach where we have many test sets.
e1 <- tsCV(mytimeseries, snaive, h=12)
e2 <- tsCV(mytimeseries, naive, h=12)
e3 <- tsCV(mytimeseries, rwf, drift=TRUE, h=12)
e4 <- tsCV(mytimeseries, meanf, h=12)
e5 <- tsCV(mytimeseries, hw,
seasonal='multiplicative', h=12)
etsfc <- function(y,h,...)
{
ets(y,...) %>% forecast(h=h)
}
e6 <- tsCV(mytimeseries, etsfc, h=12, model="MAM", damped=TRUE)
MSE <- cbind(
h=1:12,
SNaive=colMeans(tail(e1^2, -14)^2, na.rm=TRUE),
Naive=colMeans(tail(e2^2, -14)^2, na.rm=TRUE),
Drift=colMeans(tail(e3^2, -14)^2, na.rm=TRUE),
Mean = colMeans(tail(e4^2, -14)^2, na.rm=TRUE),
HW=colMeans(tail(e5^2, -14)^2, na.rm=TRUE),
ETS=colMeans(tail(e6^2, -14)^2, na.rm=TRUE))
MSE## h SNaive Naive Drift Mean HW ETS
## h=1 1 2385069 1429771 1452606 54419251 78846.56 96319.64
## h=2 2 2390801 2607311 2712352 55288993 144950.91 147594.57
## h=3 3 2396560 2502532 2558689 56154827 289634.32 264685.08
## h=4 4 2402349 3080506 3225382 57021845 473311.49 444263.15
## h=5 5 2408166 3557901 3737199 57893397 703897.35 608808.71
## h=6 6 2414011 4886609 5181524 58762022 1041599.50 726612.88
## h=7 7 2419884 4415088 4802270 59612481 1656147.78 1146700.80
## h=8 8 2425786 4294996 4717064 60489660 2433736.08 1498388.98
## h=9 9 2431714 4434867 5036448 61379960 2931203.92 1831512.04
## h=10 10 2437672 5389420 5853796 62280823 3670325.35 2438889.82
## h=11 11 2443660 5069909 5804481 63172887 4347128.21 2955566.10
## h=12 12 2449671 2449671 3026613 64094677 5191479.47 3561608.96## SNaive Naive Drift Mean HW ETS
## 1554.693 1888.516 1969.424 7692.564 1215.462 1027.048Here I have simply removed the first 14 rows to ensure we have complete forecasts for all times. It might be better to remove the first few years, to avoid the problem of having unreliable forecasts when the training data is too short.
Finally, we do the calculation allowing ets to select a different model every time.
e7 <- tsCV(mytimeseries, etsfc, h=12)
MSE <- cbind(MSE,
ETS2 = colMeans(tail(e7^2, -14)^2, na.rm=TRUE))
MSE## h SNaive Naive Drift Mean HW ETS ETS2
## h=1 1 2385069 1429771 1452606 54419251 78846.56 96319.64 96184.84
## h=2 2 2390801 2607311 2712352 55288993 144950.91 147594.57 147649.04
## h=3 3 2396560 2502532 2558689 56154827 289634.32 264685.08 264549.25
## h=4 4 2402349 3080506 3225382 57021845 473311.49 444263.15 443449.74
## h=5 5 2408166 3557901 3737199 57893397 703897.35 608808.71 607651.69
## h=6 6 2414011 4886609 5181524 58762022 1041599.50 726612.88 723696.48
## h=7 7 2419884 4415088 4802270 59612481 1656147.78 1146700.80 1144244.27
## h=8 8 2425786 4294996 4717064 60489660 2433736.08 1498388.98 1495322.60
## h=9 9 2431714 4434867 5036448 61379960 2931203.92 1831512.04 1827548.63
## h=10 10 2437672 5389420 5853796 62280823 3670325.35 2438889.82 2434087.75
## h=11 11 2443660 5069909 5804481 63172887 4347128.21 2955566.10 2944409.82
## h=12 12 2449671 2449671 3026613 64094677 5191479.47 3561608.96 3535030.89## SNaive Naive Drift Mean HW ETS ETS2
## 1554.693 1888.516 1969.424 7692.564 1215.462 1027.048 1025.460